Дз по Сетевым моделям

Подготовительные шаги

# Установим рабочую директорию
setwd("~/Documents/GitHub/RProjects_EF/Networks")

Обработка данных и создание графа

# Загрузим данные для исследования: набор ребер и вершин
edges <- read.csv("edges_1.csv", header=T) 
colnames(edges) <- c('from', 'to', 'type1', 'type2', 'weight')

nodes <- unique(edges[c('from', 'type1')])
more_nodes <- unique(edges[c('to', 'type2')])

colnames(nodes) <- c('organisation', 'type')
colnames(more_nodes) <- c('organisation', 'type')

nodes <- rbind(nodes, more_nodes)
nodes <- unique(nodes[c('organisation', 'type')])
rownames(nodes) <- NULL
edges <- edges[c('from', 'to', 'weight')]

# Создадим также набор обратных величин для весов ребер (это пригодится нам в дальнейшем). Исходные веса нам уже не понадобятся
edges$weight <- 1 / edges$weight

net <- graph_from_data_frame(d=edges, vertices=nodes, directed=F)

V(net)$color <- ifelse( V(net)$type=="Business Enterprise", "red", 
  ifelse( V(net)$type=="Private not for profit", "purple", 
    ifelse( V(net)$type=="Government", "green", 
      ifelse( V(net)$type=="Higher Education", "blue", "black") 
    ) 
  ) 
)
head(V(net)$color, 10)

 [1] "red"    "red"    "red"    "red"    "red"    "red"    "purple" "green"  "purple" "red"   

Поясним, почему в дальнейшем будем работать именно с обратными весами: Кратчайший путь и кратчайшее расстояние в данной задаче будут представлять собой не то же самое, что и в задаче с издержками перевозок или длинами расстояний между вершинами в качестве весов ребер. Отрицательные веса мы не используем, потому что алгоритм Дейкстры работает с неотрицательными весами ребер, поэтому далее мы оперируем с обратными весами, т.к., на наш взгляд, этот выбор является наиболее оптимальным для преодоления возникших трудностей.

---

Исследовательский вопрос

Исследовательский вопрос: Какие есть способы улучшения взаимодействия учреждений высшего образования Новой Зеландии (если это вообще требуется)?

Наше исследование показывает, что возможно повысить эффективность небольших университетов посредством увеличения сотрудничества с крупными, которые в свою очередь имеют тесное сотредничество с бизнесом.

Детали см. в Задании 6.

---

Анализ датасета

Описание датасета

Датасет представляет из себя множество вершин с качественными характеристиками и множество ребер с числовыми характеристиками (весами).

Вершинами являются организации Новой Зеландии, которые занимаются написанием научных статей в качестве одного из видов деятельности. Качественной характеристикой вершин является вид организации (государственная организации, высшее учебное заведение, коммерческое предприятие и частное некоммерческое предприятие). Всего в датасете 1511 вершин.

Между i-ой и j-ой вершинами (организациями) существует ребро, если за период 2010-2015 гг. на Scopus была опубликована хотя бы одна статья, хотя бы один автор которой числится в организации i, а другой автор – в организации j. Всего в датасете 4273 ребра.

Вес ребра равен количеству совместных опубликованных научных работ. Если вес ребра между i-ой и j-ой вершиной равен 3, то это означает, что за указанный период на Scopus было опубликовано 3 научные работы, у которых хотя бы один автор числится в организации i, а другой — в организации j.

---

Задание 2

При использовании нескольких пакетов для работы с графами не оказалось возможным интерпретируемо визуализировать граф полностью. Поэтому для отображения структуры нашей сети найдем 100 самых центральных по степени вершин и представим связи между ними при помощи пакета ggraph.

V(net)$degrees <- degree(net)
nodes$degrees <- V(net)$degrees

for_vis_n <- nodes[order(-nodes$degrees),][1:100, 1:2]
for_vis_e <- edges[(edges$from %in% for_vis_n$organisation & edges$to %in% for_vis_n$organisation),]

for_vis_net <- graph_from_data_frame(d=for_vis_e, vertices=for_vis_n, directed=F)
V(for_vis_net)$degrees <- degree(for_vis_net)

for_vis_net %>%
  ggraph(layout='graphopt') +
  geom_node_point(aes(color = type, size=round(degrees/20))) +
  geom_edge_link(alpha = 0.05) +
  theme_void() +
  ggtitle("Network of New Zealand organisations")

Как можно увидеть, среди 100 самых центральных вершин с количественной точки зрения преобладают организации государственного сектора, однако самыми крупными здесь являются новозеландские университеты. Первые 5 мест по числу совместных научных работ принадлежат именно им, лидером по данному показателю является Университет Окленда, самый крупный и престижный университет Новой Зеландии.

Рассчитаем базовые параметры нашего графа. Так как весами дуг в данном случае является число совместных работ, то есть “сила связи” двух вершин, будем использовать обратные веса для подсчета центральностей через кратчайшие пути. Это довольно разумно, так как от написания дополнительной общей работы прирост в силе связи двух организаций должен монотонно снижаться, как и разница в пути до соседа.

#Считаем базовые параметры графа

# Density
# The proportion of present edges from all possible ties.
p <- edge_density(net)

# связность сети
no <- igraph::components(net)$no # в пакете statnet такое же название
csize <- igraph::components(net)$csize
csize

 [1] 1463    2    2    2    3    2    2    2    2    2    2    2    2    2    2    2    2    2    2    2
[21]    2    2    3    2

# Diameter (longest geodesic distance) (если несколько компонент, то диаметр самой большой)
# Note that edge weights are used by default, unless set to NA.
diam <- diameter(net, weights=NA)
diam_w <- diameter(net)
path <- get_diameter(net, weights=NA)
path_w <- get_diameter(net)
path

+ 7/1511 vertices, named, from 83998d9:
[1] Dairy Production Systems Ltd.  Straven Road Veterinary Centre AgResearch                    
[4] Auckland Council               Auckland District Health Board New Zealand National Eye Bank 
[7] Project Litefoot Trust        

path_w

+ 9/1511 vertices, named, from 83998d9:
[1] Dairy Production Systems Ltd.     Straven Road Veterinary Centre    AgResearch                       
[4] Massey University                 University of Auckland            Victoria University of Wellington
[7] GNS Science                       ECoast Limited                    Dumpark Data Science             

# рассчитаем кратчайшие пути
# Average path length 
# The mean of the shortest distance between each pair of nodes in the network 
# (in both directions for directed graphs). 
m_dist <- mean_distance(net, weights=NA)
m_dist_w <- mean_distance(net)

Плотность нашего графа составляет 0.0037. Такая низкая плотность (вероятность образования связи) довольно естественная для вероятности написания научной работы. Наша сеть имеет 24 компонент связности: одна главная, соединяющая практически все вершины, а все остальные имеют по 2-3 вершины. Диаметр (наибольший из кратчаших путей) имеет длину 4.01, или 6, если брать единичные веса для всех дуг.

Изучим соседство самой центральной вершины — Оклендского Университета (с кем чаще всего взаимодействует наиболее престижный ВУЗ страны):

# We can also easily identify the immediate neighbors of a vertex.
# The 'neighbors' function finds all nodes one step out from the focal actor.
# To find the neighbors for multiple nodes, use 'adjacent_vertices()'.
d <- data.frame(neighbors(net, V(net)["University of Auckland"], mode="total")$type)
groups <- as.data.frame(table(d))[, 2]
pie(groups, labels=paste0(as.data.frame(table(d))[, 1], ', ',round(groups/sum(groups)*100, digits=1), '%'), 
    col=c(5, 2, 3, 4), main='Neighbors of University of Auckland')

# To find node neighborhoods going more than one step out, use function 'ego()'
# with parameter 'order' set to the number of steps out to go from the focal node(s).
neigh2 <- data.frame(ego(net, order = 2, nodes = V(net)["University of Auckland"])[[1]]$type)
groups2 <- as.data.frame(table(neigh2))[, 2]
pie(groups2, labels=paste0(as.data.frame(table(neigh2))[, 1], ', ',round(groups2/sum(groups2)*100, digits=1), '%'), 
    col=c(5, 2, 3, 4), main='Two-step neighbors of University of Auckland')

Довольно ожидаемо, что чаще всего в написании научных работ такой успешный университет сотрудничает с индустрией. Это подтвеждает наше предположение о том, успех научной работы зависит от работы с бизнесом. В ближайшем окружении данной организации также принято вести исследования с индустрией.

Сравним нашу сеть с теоритическими моделями случайного графа на основе распределения числа их соседей, диаметра, среднего кратчайшего пути и числа компонент связности. Начнем с модели случайного графа Эрдеша-Реньи.

#Генерируем модели случайного графа erdos.renyi
p <- edge_density(net)

random_model <- erdos.renyi.game(n=gorder(net), p.or.m=p, type='gnp')
par(mfrow = c(1, 2))
hist(degree(net), main = "Histogram of node degree", breaks = 1:vcount(net)-1, prob = TRUE, xlim=range(c(0, 20)))
hist(degree(random_model), main = "Histogram of node degree", breaks = 1:vcount(net)-1, prob = TRUE, xlim=range(c(0, 20)))

par(mfrow = c(1, 1))

gl <- vector('list', 100)
m <- rep(NA, 100)
d <- rep(NA, 100)
a <- rep(NA, 100)
for (i in 1:100){
  gl[[i]] <- erdos.renyi.game(n=gorder(net), p.or.m=p, type='gnp')
  m[i] <- mean_distance(gl[[i]])
  d[i] <- diameter(gl[[i]])
  a[i] <- igraph::components(gl[[i]])$no
}
mean_distance <- m
hist(mean_distance, xlim=range(c(min(m, mean_distance(net)), max(m, mean_distance(net)))))
abline(v=mean_distance(net), col='red', lty=3, lwd=2) 

diameter <- d
hist(diameter, xlim=range(c(min(d, diameter(net, weights=NA)), max(d, diameter(net, weights=NA)))))
abline(v=diameter(net, weights=NA), col='red', lty=3, lwd=2) 

number_of_components <- a
hist(number_of_components, xlim=range(c(0, 25)))
abline(v=igraph::components(net)$no, col='red', lty=3, lwd=2)

Как можно заметить, все рассматриваемые показатели нашей сети значимо отличаются от модельных, значит, структуру нашей сети нельзя описать данным образом.

Обратимся к модели “Малого мира” с регулярной решеткой. Однако в этом случае сперва нужно подобрать гиперпараметр — долю переприсоединенных дуг.

#Генерируем модели случайного графа с регулярной решеткой
n_of_nei <- mean(degree(net))
diam_net <- diameter(net, weights=NA)
#Ищем оптимальное значение p
mean_dist_net <- mean_distance(net)
model <- vector('list', 100)
m <- rep(NA, 100)
d <- rep(NA, 100)
res <- rep(NA, 100)
for (i in 1:100){
  model[[i]] <- sample_smallworld(dim = 1, size = gorder(net), nei = n_of_nei, p = i/100)
  m[i] <- mean_distance(model[[i]])
  d[i] <- diameter(model[[i]])
  res[i] <- ((m[i]-mean_dist_net)^2+(d[i]-diam_net)^2)^(1/2) 
}

# p = 0.99

random_model <- sample_smallworld(dim = 1, size = gorder(net), nei = n_of_nei, p = 0.99)
par(mfrow = c(1, 2))
hist(degree(net), main = "Histogram of node degree", breaks = 1:vcount(net)-1, prob = TRUE, xlim=range(c(0, 30)))
hist(degree(random_model), main = "Histogram of node degree", breaks = 1:vcount(net)-1, prob = TRUE, xlim=range(c(0, 30)))

par(mfrow = c(1, 1))

gl <- vector('list', 100)
m <- rep(NA, 100)
d <- rep(NA, 100)
a <- rep(NA, 100)
for (i in 1:100){
  gl[[i]] <- sample_smallworld(dim = 1, size = gorder(net), nei = n_of_nei, p = 0.99)
  m[i] <- mean_distance(gl[[i]])
  d[i] <- diameter(gl[[i]])
  a[i] <- igraph::components(gl[[i]])$no
}
mean_distance <- m
hist(mean_distance, xlim=range(c(min(m, mean_distance(net)), max(m, mean_distance(net)))))
abline(v=mean_distance(net), col='red', lty=3, lwd=2) 

diameter <- d
hist(diameter, xlim=range(c(min(d, diameter(net, weights=NA)), max(d, diameter(net, weights=NA)))))
abline(v=diameter(net, weights=NA), col='red', lty=3, lwd=2) 

number_of_components <- a
hist(number_of_components, xlim=range(c(0, 25)))
abline(v=igraph::components(net)$no, col='red', lty=3, lwd=2)

Аналогично предыдущему случаю, данная модель плохо описывает структуру нашего графа. Значимо совпадает лишь диаметр, так как по нему в том числе происходил подбор гиперпараметра.

Рассмотрим также схожесть с моделью Барабаши-Альберта, включающей предпочтительное присоединение к наиболее центральным вершинам.

# Генерируем модель Барабаши-Альберта
#Ищем оптимальное значение power
model <- vector('list', 100)
res <- rep(NA, 40)
m <- rep(NA, 100)
d <- rep(NA, 100)
for (i in c(1:40)){
  model[[i]] <- sample_pa(n = gorder(net), power = i/10, directed = FALSE)
  m[i] <- mean_distance(model[[i]])
  d[i] <- diameter(model[[i]])
  res[i] <- ((m[i]-mean_dist_net)^2+(d[i]-diam_net)^2)^(1/2) 
}

#power=2.15

random_model <- sample_pa(n = gorder(net), power = 2.15, directed = FALSE)
par(mfrow = c(1, 2))
hist(degree(net), main = "Histogram of node degree", breaks = 1:vcount(net)-1, prob = TRUE, xlim=range(c(0, 20)))
hist(degree(random_model), main = "Histogram of node degree", breaks = 1:vcount(net)-1, prob = TRUE, xlim=range(c(0, 20)))

par(mfrow = c(1, 1))

gl <- vector('list', 100)
m <- rep(NA, 100)
d <- rep(NA, 100)
a <- rep(NA, 100)
for (i in 1:100){
  gl[[i]] <- sample_pa(n = gorder(net), power = 2.15, directed = FALSE)
  m[i] <- mean_distance(gl[[i]])
  d[i] <- diameter(gl[[i]])
  a[i] <- igraph::components(gl[[i]])$no
}
mean_distance <- m
hist(mean_distance, xlim=range(c(min(m, mean_distance(net)), max(m, mean_distance(net)))))
abline(v=mean_distance(net), col='red', lty=3, lwd=2) 

diameter <- d
hist(diameter, xlim=range(c(min(d, diameter(net, weights=NA)), max(d, diameter(net, weights=NA)))))
abline(v=diameter(net, weights=NA), col='red', lty=3, lwd=2) 

number_of_components <- a
hist(number_of_components, xlim=range(c(0, 25)))
abline(v=igraph::components(net)$no, col='red', lty=3, lwd=2)

Как мы получили, данная модель тоже не подходит для описания структуры новозеландской сети: все статистически значимо отличается, кроме параметра, по которому проходил подбор гиперпараметра.

Выводы

Выводы по Заданию 2: Ни одна из теоретических моделей, известных нам, не описывает структуры нашего графа достаточно точно.

---

Задание 3

В этом задании нам предлагается посчитать различные меры центральности вершин и сделать содержательные выводы на основе полученных расчетов.

Замечание: Для анализа используем обратные веса, т.к. их использование делает осмысленным расчет центральности по близости и кратчайшему пути, но не искажает остальных мер центральности вершин.

Центральность по степени

Для начала посчитаем центральность по степени:

deg <- degree(net)

# plot(net, vertex.label=NA, layout=layout_with_fr, vertex.size=deg)

Теперь построим гистограмму распределения степеней вершин:

# гистограмма степеней вершин
hist(deg, main="Histogram of node degree", breaks=50, prob=TRUE) 

# распределение степеней вершин

mean(deg) # Средняя степень вершины по графу

[1] 5.655857

# ~ 5.66
Mode(deg) # Модальная степень вершин графа

[1] 1
attr(,"freq")
[1] 599

# 1 -> 599 раз
median(deg) # Медиана степеней вершин

[1] 2

# ~ 2
top_10_degree <- tail(sort(deg), 10) # топ-10 по степеням
top_10_degree # в топе - 7 университетов, 3 гос. учреждения

 Canterbury District Health Board                        AgResearch             University of Waikato 
                              121                               126                               136 
   Auckland District Health Board                    AUT University Victoria University of Wellington 
                              143                               145                               204 
         University of Canterbury                 Massey University               University of Otago 
                              241                               310                               406 
           University of Auckland 
                              551 

Вывод: Центральность вершин по степени не зависит от весов. Подавляющее большинство вершин имеет степень, меньшую 50. Так, медиана кол-ва соседей равняется 2. В среднем у вершины такого графа степень равна 5.66. Наиболее часто встречающаяся степень вершины — 1, которая встречается 599 раз. Если мы посмотрим 10 вершин с наибольшим числом соседей, то окажется, что 7 из них являются университетами, в числе которых крупнейший — университет Окленда (степень вершины максимальная — 551). Оставшиеся 3 организации являются государственными: 2 департамента здравоохранения и 1 исследовательская компания. Степень у всех вершин, входящих в топ-10, > 100.

Центральность по близости

Далее рассчитаем центральность по близости двумя способами: В первом варианте все веса ребер графа являются 1, в другом используются обратные веса ребер.

clos1 <- closeness(net, weights=NA) # Первый способ - без учета весов
# clos1

plot(net, vertex.label=NA, layout=layout_with_fr, vertex.size=clos1 * 50)

top_10_clos1 <- tail(sort(clos1), 10) 
# Вектор-топ 42 также будет единичным, а вот вектор-топ 43 уже нет
top_10_clos1

  Reveal Infrastructure Ltd.          Napier City Council   RATO Natural Health Clinic 
                           1                            1                            1 
Stahlton Engineered Concrete        Katipo Communications     Workwise Employment Ltd. 
                           1                            1                            1 
                  Northpower   New Zealand Stock Exchange           Real Time Genomics 
                           1                            1                            1 
    Rotary Club of Riccarton 
                           1 

bottom_10_clos1 <- head(sort(clos1), 10)
bottom_10_clos1

                                         Tauranga City Council 
                                                  0.0001571092 
                                        Project Litefoot Trust 
                                                  0.0001703868 
                                          Dumpark Data Science 
                                                  0.0001706193 
                                 Dairy Production Systems Ltd. 
                                                  0.0001720282 
                                 SciCon Scientific Consultants 
                                                  0.0001725923 
New Zealand Employment Service and Work and Income New Zealand 
                                                  0.0001773050 
                                               IWM Consultancy 
                                                  0.0001779676 
                                        PLANTwise Services Ltd 
                                                  0.0001779676 
                                   Care of Orca Research Trust 
                                                  0.0001794044 
                                    Southern Veterinary Centre 
                                                  0.0001818843 

clos1["University of Auckland"]

University of Auckland 
          0.0004175365 

clos2 <- closeness(net)
# clos2

plot(net, vertex.label=NA, layout=layout_with_fr, vertex.size=clos2 * 50)

top_10_clos2 <- tail(sort(clos2), 10) 
top_10_clos2

  Workwise Employment Ltd.                 Northpower New Zealand Stock Exchange 
                       1.0                        1.0                        1.0 
        Real Time Genomics   Rotary Club of Riccarton      Bio Soil and Crop Ltd 
                       1.0                        1.0                        1.2 
           PollenPlus Ltd.               GroPlus Ltd.           Blues Rugby Club 
                       1.2                        1.5                        2.0 
         Highlanders Rugby 
                       2.0 

bottom_10_clos2 <- head(sort(clos2), 10)
bottom_10_clos2

                                          Dumpark Data Science 
                                                  0.0002501171 
                                 Dairy Production Systems Ltd. 
                                                  0.0002502698 
New Zealand Employment Service and Work and Income New Zealand 
                                                  0.0002504864 
                                        Project Litefoot Trust 
                                                  0.0002504870 
                                         Tauranga City Council 
                                                  0.0002955979 
                                                Port Otago Ltd 
                                                  0.0003060701 
                                                 Travel Clinic 
                                                  0.0003065386 
                                         Dobbie Engineers Ltd. 
                                                  0.0003069975 
                                        NZ Strong Construction 
                                                  0.0003069975 
                                     Pulp Mill Technical Group 
                                                  0.0003069975 

clos2["University of Auckland"]

University of Auckland 
          0.0009338331 

Вывод: В обоих случаях максимальная центральность вершин будет составлять 1. В нашем графе существует 21 пара вершин, которая связана лишь между собой и более ни с кем (из рассчета для графа с единичными весами ребер). Из них лишь 5 организаций работали совместно более 1 раза (из рассчета для графа с обратными весами). Топ-10 центральных вершин по близости и топ-10 по степени уже существенно различаются. Так, центральность по близости университета Окленда составит 4.1753653^{-4} в первом случае и 4.1753653^{-4}. Т.е. round(1 / clos1["University of Auckland"]) или round(1 / clos2["University of Auckland"]) — сумма длин всех путей от университета Окленда до любой другой вершины в зависимости от метода рассчета.

Центральность по кратчайшему пути

betw1 <- betweenness(net, weights=NA)
head(betw1)

           45 South Management Ltd.                          Abacus ALS 
                             0.0000                              0.0000 
                     AbacusBio Ltd.                   ABB Power Quality 
                           228.0179                              0.0000 
ABI Rehabilitation New Zealand Ltd.              AC Research Associates 
                             0.0000                              0.0000 

top_10_betw1 <- tail(sort(betw1), 10)
top_10_betw1

   Auckland District Health Board           Plant and Food Research                        AgResearch 
                         44138.73                          55148.49                          61204.78 
            University of Waikato                    AUT University Victoria University of Wellington 
                         62949.44                          64866.37                          84858.46 
         University of Canterbury                 Massey University               University of Otago 
                        130456.08                         194334.96                         231358.36 
           University of Auckland 
                        385421.47 

bottom_10_betw1 <- head(sort(betw1), 10)
bottom_10_betw1

                45 South Management Ltd.                               Abacus ALS 
                                       0                                        0 
                       ABB Power Quality      ABI Rehabilitation New Zealand Ltd. 
                                       0                                        0 
                  AC Research Associates      Academic Colleges Group New Zealand 
                                       0                                        0 
Action on Smoking and Health New Zealand                            Active Health 
                                       0                                        0 
                         Ade Consultants                   Adidas Sports Medicine 
                                       0                                        0 

betw1["University of Auckland"]

University of Auckland 
              385421.5 

betw2 <- betweenness(net)
head(betw2)

           45 South Management Ltd.                          Abacus ALS 
                                  0                                   0 
                     AbacusBio Ltd.                   ABB Power Quality 
                                  1                                   0 
ABI Rehabilitation New Zealand Ltd.              AC Research Associates 
                                  0                                   0 

top_10_betw2 <- tail(sort(betw2), 10)
top_10_betw2

   Auckland District Health Board           Plant and Food Research                    AUT University 
                         48988.50                          57655.33                          74583.33 
            University of Waikato                        AgResearch Victoria University of Wellington 
                         81434.50                          96575.50                         139495.17 
         University of Canterbury                 Massey University               University of Otago 
                        198521.67                         323919.17                         340995.17 
           University of Auckland 
                        815924.25 

bottom_10_betw2 <- head(sort(betw2), 10)
bottom_10_betw2

                45 South Management Ltd.                               Abacus ALS 
                                       0                                        0 
                       ABB Power Quality      ABI Rehabilitation New Zealand Ltd. 
                                       0                                        0 
                  AC Research Associates      Academic Colleges Group New Zealand 
                                       0                                        0 
       Accident Compensation Corporation Action on Smoking and Health New Zealand 
                                       0                                        0 
                           Active Health                Acurity Health Group ltd. 
                                       0                                        0 

betw2["University of Auckland"]

University of Auckland 
              815924.2 

Вывод: При разных методах минимальные значения центральности у вершин будут равны 0. Это вершины, расположенные на периферии графа, или те, которые вовсе имеют всего одного соседа (например, та 21 пара вершин). Максимальные центральности вершин превышают 40 000 в топ-10. В обоих случаях университет Окленда на первом месте (3.8542147^{5} и 8.1592425^{5}). Значительное число кратчайших путей между любыми двумя вершинами проходит именно через данный университет.

Центральность по собственному значению

eig1 <- eigen_centrality(net, weights = NA)
# head(eig1)

top_10_eig1 <- tail(sort(eig1$vector), 10)
top_10_eig1

                       AgResearch    Auckland District Health Board  Canterbury District Health Board 
                        0.3760716                         0.4334652                         0.4357203 
            University of Waikato                    AUT University Victoria University of Wellington 
                        0.4373508                         0.4700463                         0.5846012 
         University of Canterbury                 Massey University               University of Otago 
                        0.5943213                         0.6779365                         0.8521192 
           University of Auckland 
                        1.0000000 

bottom_10_eig1 <- head(sort(eig1$vector), 10)
bottom_10_eig1

                                  Highlanders Rugby                                    Growth Solutionz 
                                       3.407685e-18                                        3.423996e-18 
                                           FEM Ltd. New Zealand Rock Lobster Industry Council (NZ RLIC) 
                                       3.524028e-18                                        3.563296e-18 
                                        GHD Limited                                 Napier City Council 
                                       3.577456e-18                                        3.578251e-18 
                           Rotary Club of Riccarton                              Milner Consulting Ltd. 
                                       3.585397e-18                                        3.585975e-18 
                                           NetValue                        Crown Fibre Holdings Limited 
                                       3.598880e-18                                        3.629860e-18 

eig2 <- eigen_centrality(net)
# head(eig2)

top_10_eig2 <- tail(sort(eig2$vector),10)
top_10_eig2

                             NIWA             University of Waikato                    AUT University 
                        0.1576988                         0.1897041                         0.1992335 
 Canterbury District Health Board    Auckland District Health Board Victoria University of Wellington 
                        0.2127479                         0.2346215                         0.2896466 
         University of Canterbury                 Massey University               University of Otago 
                        0.3085040                         0.4043171                         0.6707261 
           University of Auckland 
                        1.0000000 

bottom_10_eig2 <- head(sort(eig2$vector),10)
bottom_10_eig2

                            Alzheimers Northland                                  AP EnerTec Ltd. 
                                               0                                                0 
Aviation Industry Association of New Zealand Inc                            Bio Soil and Crop Ltd 
                                               0                                                0 
                                Blues Rugby Club                    Bruce Shallard and Associates 
                                               0                                                0 
                      ChangeMakers Refugee Forum                      Chatham Rock Phosphate Ltd. 
                                               0                                                0 
                      Chilton Saint James School                     Crown Fibre Holdings Limited 
                                               0                                                0 

Вывод: Минимальные значения центральностей по собственному значению также составляют 0. Наибольшую же имеет университет Окленда. Его высокий показатель говорит, что сам университет связан со многими вершинами, которые в свою очередь имеют также высокое значение центральности.

Выводы

# объединим все центральности в один датафрейм и посмотрим корреляции в двух вариантах 
df1 <- data.frame(deg, clos1, betw1, eig1$vector)
cor(df1)

                    deg       clos1       betw1 eig1.vector
deg          1.00000000 -0.03452649  0.96148582   0.8926645
clos1       -0.03452649  1.00000000 -0.01592647  -0.1091064
betw1        0.96148582 -0.01592647  1.00000000   0.7596773
eig1.vector  0.89266448 -0.10910642  0.75967732   1.0000000

df2 <- data.frame(deg, clos2, betw2, eig2$vector)
cor(df2)

                    deg       clos2       betw2 eig2.vector
deg          1.00000000 -0.03349182  0.90929438   0.8418800
clos2       -0.03349182  1.00000000 -0.01198628  -0.1332681
betw2        0.90929438 -0.01198628  1.00000000   0.7737961
eig2.vector  0.84188004 -0.13326806  0.77379607   1.0000000

Вывод: Составив корреляционные матрицы заметим, что существует высокая корреляция между следующими парами центральностей вершин: по степени - по кратчайшему расстоянию, по степени - по собств. значению и по собств. значению - по кратчайшему расстоянию

Сравним, какие организации попадали в топы по центральностям (исключая центральность по близости).

corr1 <- data.frame(
  names(top_10_degree), 
  names(top_10_betw1), 
  names(top_10_eig1)
)
corr1

corr2 <- data.frame(
  names(top_10_degree), 
  names(top_10_betw2), 
  names(top_10_eig2)
)
corr2

Выводы по Заданию 3: Сравнив топы по максимальным центральностям (за исключением топа центральностей по близости из-за малого количества соседей [порядка 1]), заметим, что по трем центральностям первую пятерку занимают одни и те же организации: университет Окленда, университет Отаго, университет Мэсси, университет Кантербери и Викторианский университет Веллингтона.

---

Задание 4

В Задании 4 нам предлагается рассчитать ассортативность и гомофилию для графа. Начнем исследование данных с ассортативности:

Ассортативность

# Возникли проблемы с вычислением ассортативности через
# assortativity_nominal(net, types=V(net)$type, directed=F) 
# Предупреждение: в результате преобразования созданы NA
# [1] NaN
assortativity_nominal(net, types=V(net)$type, directed=F)

Предупреждение: в результате преобразования созданы NA

[1] NaN

# Поэтому было принято решение создать новый столбец с закодированными строками

V(net)$num_type <- ifelse( V(net)$type=="Business Enterprise", 1, 
                          ifelse( V(net)$type=="Private not for profit", 2, 
                              ifelse( V(net)$type=="Government", 3, 
                                  ifelse( V(net)$type=="Higher Education", 4, 0) 
                              ) 
                          ) 
                       )

assort_nom <- assortativity_nominal(net, types=V(net)$num_type, directed=F) 
assort_nom # -0.05876298

[1] -0.05876298

# значение ассортативности не зависит от конкретных чисел
V(net)$num_type1 <- ifelse( V(net)$type=="Business Enterprise", 2, 
                           ifelse( V(net)$type=="Private not for profit", 7, 
                               ifelse( V(net)$type=="Government", 5, 
                                  ifelse( V(net)$type=="Higher Education", 10, 0) 
                               ) 
                           ) 
                        )

assortativity_nominal(net, types=V(net)$num_type1, directed=F) # -0.05876298

[1] -0.05876298

# считаем ассортативность по степени вершины
assort_deg <- assortativity_degree(net, directed=F) # -0.3361749
assort_deg

[1] -0.3361749

Обобщение результатов: Ассортативность в графе по типу вершин примерно равна -0.059. Это означает, что мы не можем утверждать наличие предпочтительного присоединения узлов сети к своему типу (среди узлов в общем, т.е. среди узлов всех типов).

Ассортативность по степени вершины составляет примерно -0.336, что говорит о предпочтительном присоединении узлов к вершинам с относительно небольшим количеством соседей.

Для выявления более точного эффекта предпочтительного присоединения именно среди университетов (что интересует нас в рамках исследовательского вопроса) рассчитаем показатели гомофилии для университетов и других типов институтов.

Нами было принято решение рассчитать dyadicity (ее посчитаем только для одного случая, т.к. она не будет меняться при различных разбиениях) и heterophilicity учреждений высшего образования для следующих разделений на подгруппы: - Higher Education vs все остальные - Higher Education vs Business Enterprise - Higher Education vs Government.

Замечание: Влиянием Private not for profit можно пренебречь.

Higher Education vs все остальные

p <- edge_density(net)
p # 0.003745601

[1] 0.003745601

n_c <- sum( V(net)$type == 'Higher Education' )
n_c # 54

[1] 54

exp_diad_uni <- n_c * (n_c - 1) / 2 * p
exp_diad_uni # 5.359955

[1] 5.359955

unis <- V(net)[which( V(net)$type == 'Higher Education' )]
unis <- names(unis)
# unis

m_c <- 0

for (row in 1:nrow(edges)) {
  from_obj <- edges[row, 'from']
  to_obj <- edges[row, 'to']
  
  if ( (from_obj %in% unis) && (to_obj %in% unis) ) {
    m_c <- m_c + 1
  }
}

m_c # 172

[1] 172

diadicity_uni <- m_c / exp_diad_uni 
diadicity_uni # 32.08982

[1] 32.08982

# 2-й способ посчитать dyadicity
uni_net <- induced_subgraph( 
  net, 
  v=V(net)[which( V(net)$type == 'Higher Education' )] 
)
# uni_net
# V(uni_net)$type # проверяем, что все вершины в новом графе вузы

edge_density(net) # 0.003745601

[1] 0.003745601

edge_density(uni_net) # 0.1201957

[1] 0.1201957

edge_density(uni_net) / edge_density(net) # 32.08982

[1] 32.08982

n_c # 54

[1] 54

p # 0.003745601

[1] 0.003745601

n_nc <- sum( V(net)$type != 'Higher Education' )
n_nc # 1457

[1] 1457

exp_hetero_uni <- n_nc * n_c * p
exp_hetero_uni # 294.6964

[1] 294.6964

m_mixed <- 0

for (row in 1:nrow(edges)) {
  from_obj <- edges[row, 'from']
  to_obj <- edges[row, 'to']
  
  if ( (from_obj %in% unis) && !(to_obj %in% unis) ) {
    m_mixed <- m_mixed + 1
  }
  if ( !(from_obj %in% unis) && (to_obj %in% unis) ) {
    m_mixed <- m_mixed + 1
  }
}

m_mixed # 2051

[1] 2051

heterophilicity_uni <- m_mixed / exp_hetero_uni
heterophilicity_uni # 6.959706

[1] 6.959706

Higher Education vs Business Enterprise

# uni_biz_net <- induced_subgraph( 
#   net, 
#   v=V(net)[which( 
#     (V(net)$type == 'Higher Education') || (V(net)$type == 'Business Enterprise') 
#       )] 
# )
uni_biz_net <- induced_subgraph( 
  net, 
  v=V(net)[which( V(net)$type %in% c('Higher Education', 'Business Enterprise') )] 
)
# uni_biz_net
# V(uni_biz_net)$type

uni_wout_biz_net <- induced_subgraph( 
  uni_biz_net, 
  v=V(uni_biz_net)[which( V(uni_biz_net)$type == 'Higher Education' )] 
)
# uni_wout_biz_net
# V(uni_wout_biz_net)$type


p_uni_vs_biz <- edge_density(uni_biz_net)
p_uni_vs_biz # 0.003036165

[1] 0.003036165

edge_density(uni_wout_biz_net) # 0.1201957

[1] 0.1201957

dyadicity_uni_vs_biz <- edge_density(uni_wout_biz_net) / p_uni_vs_biz 
dyadicity_uni_vs_biz # 39.58799

[1] 39.58799

n_c # 54

[1] 54

p_uni_vs_biz # 0.003036165

[1] 0.003036165

n_nc_biz <- sum( V(net)$type == 'Business Enterprise' )
n_nc_biz # 1010

[1] 1010

exp_hetero_uni_vs_biz <- n_nc_biz * n_c * p_uni_vs_biz
exp_hetero_uni_vs_biz # 165.5925

[1] 165.5925

biz <- V(net)[which( V(net)$type == 'Business Enterprise' )]
biz <- names(biz)
# biz

m_mixed_vs_biz <- 0

for (row in 1:nrow(edges)) {
  from_obj <- edges[row, 'from']
  to_obj <- edges[row, 'to']
  
  if ( (from_obj %in% unis) && (to_obj %in% biz) ) {
    m_mixed_vs_biz <- m_mixed_vs_biz + 1
  }
  if ( (from_obj %in% biz) && (to_obj %in% unis) ) {
    m_mixed_vs_biz <- m_mixed_vs_biz + 1
  }
}

m_mixed_vs_biz # 1139

[1] 1139

heterophilicity_uni_vs_biz <- m_mixed_vs_biz / exp_hetero_uni_vs_biz
heterophilicity_uni_vs_biz # 6.878333

[1] 6.878333

Higher Education vs Government

uni_gov_net <- induced_subgraph( 
  net, 
  v=V(net)[which( V(net)$type %in% c('Higher Education', 'Government') )] 
)
# V(uni_gov_net)$type

uni_wout_gov_net <- induced_subgraph( 
  uni_gov_net, 
  v=V(uni_gov_net)[which( V(uni_gov_net)$type == 'Higher Education' )] 
)
# V(uni_wout_gov_net)$type


p_uni_vs_gov <- edge_density(uni_gov_net) 
p_uni_vs_gov # 0.04231381

[1] 0.04231381

edge_density(uni_wout_gov_net) # 0.1201957

[1] 0.1201957

dyadicity_uni_vs_gov <- edge_density(uni_wout_gov_net) / p_uni_vs_gov 
dyadicity_uni_vs_gov # 2.840578

[1] 2.840578

n_c # 54

[1] 54

p_uni_vs_gov # 0.04231381

[1] 0.04231381

n_nc_gov <- sum( V(net)$type == 'Business Enterprise' )
n_nc_gov # 1010

[1] 1010

exp_hetero_uni_vs_gov <- n_nc_gov * n_c * p_uni_vs_gov
exp_hetero_uni_vs_gov # 2307.795

[1] 2307.795

gov <- V(net)[which( V(net)$type == 'Business Enterprise' )]
gov <- names(gov)
# gov

m_mixed_vs_gov <- 0

for (row in 1:nrow(edges)) {
  from_obj <- edges[row, 'from']
  to_obj <- edges[row, 'to']
  
  if ( (from_obj %in% unis) && (to_obj %in% gov) ) {
    m_mixed_vs_gov <- m_mixed_vs_gov + 1
  }
  if ( (from_obj %in% gov) && (to_obj %in% unis) ) {
    m_mixed_vs_gov <- m_mixed_vs_gov + 1
  }
}

m_mixed_vs_gov # 0.4935447

[1] 1139

heterophilicity_uni_vs_gov <- m_mixed_vs_gov / exp_hetero_uni_vs_gov
heterophilicity_uni_vs_gov # 0.4935447

[1] 0.4935447

Соберем полученные значения dyadicity и heterophilicity вместе: - Higher Education vs все остальные: - dyadicity: 32.0898241 - heterophilicity: 6.9597056 - Higher Education vs Business Enterprise: - dyadicity: 39.5879866 - heterophilicity: 6.8783329 - Higher Education vs Government: - dyadicity: 2.840578 - heterophilicity: 0.4935447.

Мы видим, что dyadicity университетов при разделении на университеты и всех остальных (а также при разделении на университеты и бизнес) очень высокая. При этом heterophilicity также доволльно сильно превышает 1.

В чем может быть объяснение этой особенности данных? Как утверждается в статье “Functional dyadicity and heterophilicity of gene-gene interactions in statistical epistasis networks”, описанный выше факт может быть связан с тем, что вершины со значением метки 1 имеют большую центральность \(\Rightarrow\) хорошо связаны друг с другом и с вершинами другого типа.

В Задании 6 мы будем проверять гипотезу касательно высоких степеней вершин типа Higher Education среди крупных университетов.

Выводы

Выводы по Заданию 4: Мы получили, что несмотря на плотные связи учреждений высшего образования между собой, по крайней мере некоторые из них имеют также прочные связи с другими типами организаций (в частности с бизнесом).

---

Задание 5

Проведем процедуру кластеризации вершин при помощи Edge-betweenness и Fast-greedy методов и сравним полученные результаты с существующей классификацией.

Т.к. считать кратчайший путь для исходных весов ребер является не очень осмысленным, используем обратные для исходных веса ребер.

Замечание: Для Fast-greedy не важно, какие данные использовать: исходные или инвертированные, т.к. этот метод не учитывает веса ребер (только их количество).

# проводим кластеризацию при помощи метода Edge-betweenness
ceb <- cluster_edge_betweenness(net, weights=E(net)$weight)

Предупреждение: At core/community/edge_betweenness.c:493 : Membership vector will be selected based on the highest modularity score.

# смотрим на полученные в результате кластеризации объекты
mbship_ceb <- membership(ceb) 
# mbship_ceb

s_ceb <- sizes(ceb)
s_ceb

Community sizes
  1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26 
 45 213 233 154 229 190  75  46 164  48  32   2   2   2   3  15   2   2   8   2   2   4   2   2   2   2 
 27  28  29  30  31  32  33  34  35  36  37  38 
  7   2   2   2   2   2   2   2   2   2   3   2 

ms_ceb <- modularity(ceb)
ms_ceb # 0.527262

[1] 0.527262

# ceb

# > s_ceb <- sizes(ceb)
# > s_ceb
# Community sizes
#   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20 
#  45 213 233 154 229 190  75  46 164  48  32   2   2   2   3  15   2   2   8   2 
#  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38 
#   2   4   2   2   2   2   7   2   2   2   2   2   2   2   2   2   3   2 

# > ms_ceb <- modularity(ceb)
# > ms_ceb
# [1] 0.527262

Визуализация кластеризации методом Edge-betweenness:

plot(
  ceb, net, 
  vertex.label=NA, 
  vertex.color=V(net)$color, 
  layout=layout_nicely
)

plot(
  ceb, net, 
  vertex.label=NA, 
  vertex.color=V(net)$color, 
  layout=layout_with_lgl
)

Сравним результаты кластеризации Edge-betweenness и имеющуюся классификацию вершин:

compare(ceb, V(net)$type, method='nmi') # 

[1] 0.03399001

# здесь и далее мы используем нормализованное значение функции compare() для облегчения интерпретации возвращаемого результата

Таким образом, мы видим, что имеющаяся классификация вершин совсем не похожа на то, что мы получили при кластеризации. По всей видимости, “разделение по бутылочному горлышку” является не очень хорошим способом для кластеризации вершин в нашем датасете.

# проводим кластеризацию при помощи метода Fast-greedy
cfg <- cluster_fast_greedy(net, weights=E(net)$weight) 

# смотрим на полученные в результате кластеризации объекты
mbship_cfg <- membership(cfg) 
# mbship_cfg

s_cfg <- sizes(cfg)
s_cfg

Community sizes
  1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26 
219 144 212 132 141  94 189 120  84  82  28   3   2  15   3   3   2   2   2   2   2   2   2   2   2   2 
 27  28  29  30  31  32  33  34  35  36 
  2   2   2   2   2   2   2   2   2   2 

ms_cfg <- modularity(cfg)
ms_cfg # 0.5310055

[1] 0.5310055

# cfg

# > s_cfg <- sizes(cfg)
# > s_cfg
# Community sizes
#   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26 
# 219 144 212 132 141  94 189 120  84  82  28   3   2  15   3   3   2   2   2   2   2   2   2   2   2   2 
#  27  28  29  30  31  32  33  34  35  36 
#   2   2   2   2   2   2   2   2   2   2 

# > ms_cfg <- modularity(cfg)
# > ms_cfg
# [1] 0.5310055

Визуализация кластеризации методом Fast-greedy:

plot(
  cfg, net, 
  vertex.label=NA, 
  vertex.color=V(net)$color, 
  layout=layout_nicely
)

plot(
  cfg, net, 
  vertex.label=NA, 
  vertex.color=V(net)$color, 
  layout=layout_with_lgl
)

compare(cfg, V(net)$type, method='nmi') # 0.03386511

[1] 0.03386511

Таким образом, мы видим, что так же, как и предыдущий способ кластеризации, Fast-greedy не очень хорошо соответствует реальной классификации рассматриваемых институтов.

Выводы

Вывод по Заданию 5: Известные нам способы кластеризации (Edge-betweenness и Fast-greedy) не дали результата, похожего на имеющуюся в датасете классификацию.

---

Задание 6

Изначально был поставлен исследовательский вопрос, насколько активно при написании научных работ университеты Новой Зеландии взаимодействуют с другими типами организаций, нужно ли им как-то в этом способстовать извне (например, создавать общие лаботарии с индустриальными корпорациями и государственными органами). Помимо расчета значений гомофилии и различных видов центральностей, для ответа на заданный вопрос нам может помочь модель предсказания класса вершин Relational Neighbors Classifier (RNC).

Разделим все вершины на 2 класса: “1” - университет, “0” - иначе. Отберем из каждого класса по 10 вершин и положим вероятности их принадлежности первому классу соответственно 1 и 0. Для всех остальных вершин (unknown) положим эту вероятность равной 0.5. При каждой итерации будем обновлять вероятности каждой вершины исходя из вероятностей её соседей. В итоге, после 10 таких итераций получим для каждой вершины вероятность её принадлежности классу “1”, и отберем те вершины, чья вероятность больше порога 0.5 в качестве принадлежащих классу “1”.

Для избежания случайности повторим этот алгоритм 100 раз и построим гистограмму распределения числа предсказанных университетов, а также гистограмму распределения средних вероятностей принадлежности группе университетов.

res <- rep(NA, 20)
pr <- rep(NA, 20)

for (j in 1:20){ 
  network <- net
  V(network)[which( V(network)$type == 'Higher Education' )]$type = 1
  V(network)[which(V(network)$type != 1)]$type = 0
  
  
  # вероятность вершины перейти в класс churn = 1 равна доле соседей churn = 1 от общего числа соседей
  # используем матричные вычисления
  random_uni <- sample( V(network)[ which(V(network)$type == 1) ], 10 )
  random_not_uni <- sample( V(network)[ which(V(network)$type == 0) ], 10 )
  
  
  for (i in 1:1511){
    if (!(V(network)[i] %in% random_uni | (V(network)[i] %in% random_not_uni))){
      V(network)$type[i] = 0.5
    }
  }
  
  UniProb <-  as.numeric(V(network)$type)
  threshold <- 0.5
  A <- as_adjacency_matrix(network)
  A <- as.matrix(A)
  Neighbors <- rowSums(A) # общее число соседей вершины
  
  UniProb_iter <- UniProb
  for(i in 1:10){
    UniProb_iter <- as.vector((A %*% UniProb_iter) / Neighbors)
  }
  res[j] <- length(which(UniProb_iter > threshold))
  pr[j] <- mean(UniProb_iter)
}

hist(res, main='Число предсказанных университетов')

hist(pr, main='Средняя вероятность принадлежности классу 1')

Как можно увидеть, почти во всех случаях абсолютное число вершин предсказанj в качестве университетов. Получаем, что какие бы 10 университетов мы ни взяли, их класс дойдет до большинства вершин быстрее, чем противоположный класс 10 случайных вершин. Это распространяется и на случай, когда были выбраны наименее центральные университеты. Если бы университеты были плохо связаны с организациями другого типа, то их класс намного слабее распространялся бы на вершины другого типа, а это не так.

Значит, по крайней мере какие-то университеты (а именно наиболее центральные) довольно активно кооперируются с другими организациями (как, например, это делает Оклендский университет). Поэтому для таких университетов нет необходимости строить совместные с бизнесом и государством научные центры.

Вспомним про высокий показатель dyadicity у ВУЗов: выходит, что он получается таким большим из-за высокой гомофилии более мелких учебных заведений. Тогда можно предположить, что именно они плохо взаимодействуют с бизнесом и гос. органами.

---

Заключение

Как тогда можно им помочь? При моделировании их класс распространялся по сети через их соседей — крупные университеты. В таком случае, более разумно строить научные центры именно соединяющие мелкие университеты с крупными. Данным образом, казалось бы, усиливая гомофилию, мы, наоборот, будем способстовать развитию научных отношений университетов с индустриальными партнерами.

Замечание

Здесь мы предполагаем эталонность крупных университетов (например, Оклендского) с точки зрения эффективности взаимодействия, а также качества и количества научных работ, что продемонстрировано, например, в этой статье.